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(57) Abstract: In a solid tumor or other cancerous 
tissue in a human or animal patient, specific objects or 
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and its progress. In a three-dimensional image of the 
region of interest (Fig .1 element 102), the biomarkers 
are identified and quantified (Fig. elements 104). 
Multiple three-dimensional images can be taken over 
time, in which the biomarkers can be tracked over time 
(Fig. element 108). Statistical segmentation techniques 
are used to identify the biomarker in a first image to 
carry the identification over to the remaining images. 
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SYSTEM AND METHOD FOR QUANTITATIVE ASSESSMENT OF CANCERS 
AND THEIR CHANGE OVER TIME 
Reference to Related Applications 

The present application claims the benefit of U.S. Provisional Application No. 

5 60/322,427, filed September 17, 2001, whose disclosure is hereby incorporated by reference 
in its entirety into the present disclosure. 
Field of the Invention 

The present invention is directed to a system and method for quantifying cancers and 
their change over time and is more particularly directed to such a system and method which 

10 use biomarlcers related to cancers, or oncomarkers. 
Description of Related Art 

Malignant tumors, including cancers of the lungs, abdominal organs, bones, and 
central nervous system, afflict a significant percent of the population. In assessing those 
conditions, and in tracking their change over time, including improvements due to new 

15 therapies, it is necessary to have quantitative information. Manually obtained and imprecise 
measures of tumor growth, traditionally assessed through manual tracings or by caliper 
measurements of an image, have been used in the past. Such measures lack sensitivity and 
are typically useful only for gross characterization of tumor behavior. Examples of 
measurements that are taken from MRI or CT examinations of cancer patients include: lesion 

20 volume, lesion surface area within one slice, major and minor axes within one slice, and the 
cross product of major and minor axes within one slice. 

Some references for the prior work include: Therasse, P., et al. "New Guidelines to 
Evaluate the Response to Treatment in Solid Tumors," Journal of National Cancer Institute, 
Feb. 2000(92)3: 205-216. That paper describes the standard (RECIST) for unidimensional 

25 tumor measurement. 
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Also, for an example of the awkwardness of the conventional mouse-driven manual 
outlining of lesions, see: Barseghian, T. "Uterine Fibroid Embolization Offers Alternative to 
Surgery ," Diagnostic Imaging, Sept. 1997, 11-12. 

Other references include: 
5 Pieterman, R. et al. "Preoperative Staging of Non-Small-Cell Lung Cancer with 

Positron-Emission Tomography," New England Journal of Medicine, 2000 Jul. 27, 343(4) 
290-2. 

Yang, W., et al., "Comparison of Dynamic Helical CT and Dynamic MR Imaging in 
the Evaluation of Pelvic Lymph Nodes in Cervical Carcinoma, " American Journal of 
10 Roentgenology, 2000 Sep; 175(3) 759 - 766. 

Lilleby, W., et al. "Computed Tomography/Magnetic Resonance Based Volume 
Changes of the Primary Tumour in Patients with Prostate Cancer with or without Androgen 
Deprivation," Radiotherapy and Oncology, 2000 Nov; 57(2): 195-200. 

Ward, R,, et al. "Phase I Clinical Trial of the Chimeric Monoclonal Antibody (C30.6) 
15 in Patients with Metastatic Colorectal Cancer," Clinical Cancer Research, 2000 Dec; 6(12): 
4674-4683. 

Hermans, R., et al. "The Relation of CT-Determined Tumor Parameters and Local and 
Regional Outcome of Tonsillar Cancer after Definitive Radiation Treatment," International 
Journal of Radiation Oncology Biology-Physics. 2001 May 1; 50(1): 37-45. 
20 Stbkkel, M., et al. "Staging of Lymph Nodes with FDG Dual-Headed PET in Patients 

with Non-Small-Cell Lung Cancer," Nuclear Medicine Communications, 1 999 Nov; 
20(11):1001-1007. 

Sahani, D., et al. "Quantitative Measurements of Medical Images for Pharmaceutical 
Clinical Trials: Comparison Between On and Off-Site Assessments," American Journal of 
25 Roentgenology, 2000 Apr; 174(4): 1159-1162. 
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Couteau, C, et al., "A Phase II Study of Docetaxel in Patients with Metastatic 
Squamous Cell Carcinoma of the Head and Neck" British Journal of Cancer, 1999 Oct; 
81(3):457-462. 

Padhani, A., et al. "Dynamic Contrast Enhanced MRI of Prostate Cancer: Correlation 
5 with Morphology and Tumour Stage, Histologic Grade and PSA," Clinical Radiology, 2000 
Feb; 55(2): 99-109. 

Yankelevitz, D., et al. "Small Pulomonary Nodules: Volumetrically Determined 
Growth Rates Based on CT Evaluation," Radiology, 2000 Oct; 217: 251-256. 

Those measurements require manual or semi-manual systems that require a user to 
10 identify the structure of interest and to trace boundaries or areas, or to initialize an active 
contour. 

The prior art is capable of assessing gross changes over time. However, the 
conventional measurements are not well suited to assessing and quantifying subtle changes in 
lesion size, and are incapable of describing complex topology or shape in an accurate manner 
15 or of addressing finer details of tumor biology. Furthermore, manual and semi-manual 
measurements from raw images suffer from a high inter-observer and intra-observer 
variability. Also, manual and semi-manual measurements tend to produce ragged and 
irregular boundaries in 3D when the tracings are based on a sequence of 2D images. 
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Summary of the Invention 

It will be apparent from the above that a need exists in the art to identify features of 
tumors such as their boundaries and sub-components. It is therefore a primary object of the 
invention to provide a more accurate quantification of solid tumors and other cancerous 
5 tissues. It is another object of the invention to provide a more accurate quantification of 
changes in time of those tissues. It is a further object of the invention to address the needs 
noted above. 

To achieve the above and other objects, the present invention is directed to a 
technique for identifying characteristics of cancerous tissue, such as tumor margins, and 

10 identifying specific sub-components such as necrotic core, viable perimeter, and development 
of tumor vasculature (angiogenesis), which are sensitive indicators of disease progress or 
response to therapy. The topological, morphological, radiological, and pharmacokinetic 
characteristics of tumors and their sub-structures are called biomarkers, and specific 
measurements of the biomarkers serve as the quantitative assessment of disease progress. 

15 Biomarkers specific to tumors are also called oncomarkers. 

The inventors have discovered that the following new biomarkers are sensitive 
indicators of the progress of diseases characterized by solid tumors and other cancerous 
tissues in humans and in animals: 

• tumor surface area; 

20 • tumor compactness (surface-to-volume ratio); 

• tumor surface curvature; 

• tumor surface roughness; 

• necrotic core volume; 

• necrotic core compactness; 
25 • necrotic core shape; 
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viable periphery volume; 
volume of tumor vasculature; 
change in tumor vasculature over time; 
tumor shape, as defined through spherical harmonic analysis; 
morphological surface characteristics; 
lesion characteristics; 
tumor characteristics; 
tumor peripheral characteristics; 
tumor core characteristics; 
10 • bone metastases characteristics; 
ascites characteristics; 
pleural fluid characteristics; 
vessel structure characteristics; 
neovasculature characteristics; 
15 • polyp characteristics; 

nodule characteristics; 
angiogenisis characteristics; 
tumor length; 
tumor width; and 
20 • tumor 3d volume. 

A preferred method for extracting the biomarkers is with statistical based reasoning as 
defined in Parker et al (US Patent 6,169,817), whose disclosure is hereby incorporated by 
reference in its entirety into the present disclosure. A preferred method for quantifying shape 
and topology is with the morphological and topological formulas as defined by the following 
25 references: 

5 
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Curvature Analysis: Peet, F.G., Sahota, T.S. "Surface Curvature as a Measure of 
Image Texture" IEEE Transactions on Pattern Analysis and Machine Intelligence 1985 Vol 
PAMI-7 G:734-738; 

Struik, D J., Lectures on Classical Differential Geometry, 2nd ed., Dover, 1988. 
5 Shape and Topological Descriptors: Duda, R.O, Hart, P.E., Pattern Classification and 

Scene Analysis, Wiley & Sons, 1973. 

Jain, A.K, Fundamentals of Digital Image Processing, Prentice Hall, 1989. 

Spherical Harmonics: Matheny, A., Goldgof, D. "The Use of Three and Four 
Dimensional Surface Harmonics for Nonrigid Shape Recovery and Representation," IEEE 
10 Transactions on Pattern Analysis and Machine Intelligence 1995, 17: 967-981 ; 

Chen, C.W, Huang, T.S., Arrot, M. "Modeling, Analysis, and Visualization of Left 
Ventricle Shape and Motion by Hierarchical Decomposition," IEEE Transactions on Pattern 
Analysis and Machine Intelligence 1994, 342-356. 

Those morphological and topological measurements have not in the past been applied 
15 to onco-biomarkers. 

The quantitative assessment of the new biomarkers listed above provides an objective 
measurement of the state of progression of diseases characterized by solid tumors. It is also 
very useful to obtain accurate measurements of those biomarkers over time, particularly to 
judge the degree of response to a new therapy. Manual and semi-manual assessments of 
20 conventional biomarkers (such as major axis length or cross-sectional area) have a high 
inherent variability, so that as successive scans are traced, the variability can hide subtle 
trends. That means that only gross changes, sometimes over very long time periods, can be 
verified using conventional methods. The inventors have discovered that extracting the 
biomarker using statistical tests and treating the biomarker over time as a 4D object, with an 
25 automatic passing of boundaries from one time interval to the next, can provide a highly 
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accurate and reproducible segmentation from which trends over time can be detected. That 
preferred approach is defined in the above-cited patent to Parker et al. Thus, the combination 
of selected biomarkers that themselves capture subtle pathologies with a 4D approach to 
increase accuracy and reliability over time, creates sensitivity that has not been previously 
5 obtainable. 

The quantitative measure of the tumor can be one or more of tumor shape, tumor 
surface morphology, tumor surface curvature and tumor surface roughness. 
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Brief Description of the Drawings 

A preferred embodiment of the present invention will be set forth in detail with 
reference to the drawings, in which: 

Fig. 1 shows a flow chart of an overview of the process of the preferred embodiment; 

Fig. 2 shows a flow chart of a segmentation process used in the process of Fig. 1; 

Fig. 3 shows a process of tracking a segmented image in multiple images taken over 
time; and 

Fig. 4 shows a block diagram of a system on which the process of Figs. 1-3 can be 
implemented. 
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Detailed Description of the Preferred Embodiment 

A preferred embodiment of the present invention will now be set forth with reference 
to the drawings. 

Fig. 1 shows an overview of the process of identifying biomarkers and their trends 
5 over time. In step 102, a three-dimensional image of the region of interest is taken. In step 
104, at least one biomarker is identified in the image; the technique for doing so will be 
explained with reference to Fig. 2. Also in step 104, at least one quantitative measurement is 
made of the biomarker. In step 106, multiple three-dimensional images of the same region of 
the region of interest are taken over time. In some cases, step 106 can be completed before 

10 step 104; the order of the two steps is a matter of convenience. In step 108, the same 
biomarker or biomarkers and their quantitative measurements are identified in the images 
taken over time; the technique for doing so will be explained with reference to Fig. 3. The 
identification of the biomarkers in the multiple image allows the development in step 1 10 of a 
model of the region of interest in four dimensions, namely, three dimensions of space and one 

15 of time. From that model, the development of the biomarker or biomarkers can be tracked 
over time in step 112. 

The preferred method for extracting the biomarkers is with statistical based reasoning 
as defined in Parker et al (US Patent 6,169,817), whose disclosure is hereby incorporated by 
reference in its entirety into the present disclosure. From raw image data obtained through 

20 magnetic resonance imaging or the like, an object is reconstructed and visualized in four 
dimensions (both space and time) by first dividing the first image in the sequence of images 
into regions through statistical estimation of the mean value and variance of the image data 
and joining of picture elements (voxels) that are sufficiently similar and then extrapolating 
the regions to the remainder of the images by using known motion characteristics of 

25 components of the image (e.g., spring constants of muscles and tendons) to estimate the rigid 
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and deformational motion of each region from image to image. The object and its regions 
can be rendered and interacted with in a four-dimensional (4D) virtual reality environment, 
the four dimensions being three spatial dimensions and time. 

The segmentation will be explained with reference to Fig. 2. First, at step 201, the 

5 images in the sequence are taken, as by an MRI. Raw image data are thus obtained. Then, at 
step 203, the raw data of the first image in the sequence are input into a computing device. 
Next, for each voxel, the local mean value and region variance of the image data are 
estimated at step 205. The connectivity among the voxels is estimated at step 207 by a 
comparison of the mean values and variances estimated at step 205 to form regions. Once the 

10 connectivity is estimated, it is determined which regions need to be split, and those regions 
are split, at step 209. The accuracy of those regions can be improved still more through the 
segmentation relaxation of step 211. Then, it is determined which regions need to be merged, 
and those regions are merged, at step 213. Again, segmentation relaxation is performed, at 
step 215. Thus, the raw image data are converted into a segmented image, which is the end 

15 result at step 217. Further details of any of those processes can be found in the above-cited 
Parker et al patent. 

The creation of a 4D model (in three dimensions of space and one of time) will be 
described with reference to Fig. 3. A motion tracking and estimation algorithm provides the 
information needed to pass the segmented image from one frame to another once the first 
20 image in the sequence and the completely segmented image derived therefrom as described 
above have been input at step 301. The presence of both the rigid and non-rigid components 
should ideally be taken into account in the estimation of the 3D motion. According to the 
present invention, the motion vector of each voxel is estimated after the registration of 
selected feature points in the image. 
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To take into consideration the movement of the many structures present in the region 
of interest, the approach of the present invention takes into account the local deformations of 
soft tissues by using a priori knowledge of the material properties of the different structures 
found in the image segmentation. Such knowledge is input in an appropriate database form at 
5 step 303. Also, different strategies can be applied to the motion of the rigid structures and to 
that of the soft tissues. Once the selected points have been registered, the motion vector of 
every voxel in the image is computed by interpolating the motion vectors of the selected 
points. Once the motion vector of each voxel has been estimated, the segmentation of the 
next image in the sequence is just the propagation of the segmentation of the former image. 

10 That technique is repeated until every image in the sequence has been analyzed. The 
definition of time and the order of a sequence can be reversed for convenience in the analysis. 

Finite-element models (FEM) are known for the analysis of images and for time- 
evolution analysis. The present invention follows a similar approach and recovers the point 
correspondence by minimizing the total energy of a mesh of masses and springs that models 

15 the physical properties of the anatomy. In the present invention, the mesh is not constrained 
by a single structure in the image, but instead is free to model the whole volumetric image, in 
which topological properties are supplied by the first segmented image and the physical 
properties are supplied by the a priori properties and the first segmented image. The motion 
estimation approach is an FEM-based point correspondence recovery algorithm between two 

20 consecutive images in the sequence. Each node in the mesh is an automatically selected 
feature point of the image sought to be tracked, and the spring stiffness is computed from the 
first segmented image and a priori knowledge of the human anatomy and typical 
biomechanical properties for the tissues in the region of interest. 

Many deformable models assume that a vector force field that drives spring-attached 

25 point masses can be extracted from the image. Most such models use that approach to build 
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semi-automatic feature extraction algorithms. The present invention employs a similar 
approach and assumes that the image sampled at t = n is a set of three dynamic scalar fields: 
<D(x,0=te*(x), |Vg„(x)|,V 2 ^(x)}, 

namely, the gray-scale image value, the magnitude of the gradient of the image value, and the 
5 Laplacian of the image value. Accordingly, a change in <D(x, t) causes a quadratic change in 
the scalar field energy Uq(x) qc (AO(x)) 2 . Furthermore, the structures underlying the image 
are assumed to be modeled as a mesh of spring-attached point masses in a state of 
equilibrium with those scalar fields. Although equilibrium assumes that there is an external 
force field, the shape of the force field is not important. The distribution of the point masses 
10 is assumed to change in time, and the total energy change in a time period At after time t = n 
is given by 
AU tt (Ax) = 

I Mg m (x)-g M {x + Ax))) 2 +(A|V^W|-|Vg w+1 (^ + Ax)|)) 2 + 
(y(V a ft(*) + V 2 g w+1 (* + Ax))) 2 +±?jAX T KAX} 

where a, p, and y are weights for the contribution of every individual field change, r\ weighs 
the gain in the strain energy, K is the FEM stifihess matrix, and AX is the FEM node 
15 displacement matrix. Analysis of that equation shows that any change in the image fields or 
in the mesh point distribution increases the system total energy. Therefore, the point 
correspondence from g„ to g n +\ is given by the mesh configuration whose total energy 
variation is a minimum. Accordingly, the point correspondence is given by 

X = X + hX 
20 where 

AA^min^ AU n (AX). 

In that notation, min^ q is the value of p that minimizes q. 

12 



WO 03/025837 PCT7US02/29005 

While the equations set forth above could conceivably be used to estimate the motion 
(point correspondence) of every voxel in the image, the number of voxels, which is typically 
over one million, and the complex nature of the equations make global minimization difficult. 
To simplify the problem, a coarse FEM mesh is constructed with selected points from the 
5 image at step 305. The energy minimization gives the point correspondence of the selected 
points. 

The selection of such points is not trivial. First, for practical purposes, the number of 
points has to be very small, typically = 10 4 ; care must be taken that the selected points 
describe the whole image motion. Second, region boundaries are important features because 

10 boundary tracking is enough for accurate region motion description. Third, at region 
boundaries, the magnitude of the gradient is high, and the Laplacian is at a zero crossing 
point, making region boundaries easy features to track. Accordingly, segmented boundary 
points are selected in the construction of the FEM. 

Although the boundary points represent a small subset of the image points, there are 

15 still too many boundary points for practical purposes. In order to reduce the number of 
points, constrained random sampling of the boundary points is used for the point extraction 
step. The constraint consists of avoiding the selection of a point too close to the points 
already selected. That constraint allows a more uniform selection of the points across the 
boundaries. Finally, to reduce the motion estimation error at points internal to each region, a 

20 few more points of the image are randomly selected using the same distance constraint. 
Experimental results show that between 5,000 and 10,000 points are enough to estimate and 
describe the motion of a typical volumetric image of 256x256x34 voxels. Of the selected 
points, 75% are arbitrarily chosen as boundary points, while the remaining 25% are interior 
points. Of course, other percentages can be used where appropriate, 
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Once a set of points to track is selected, the next step is to construct an FEM mesh for 
those points at step 307. The mesh constrains the kind of motion allowed by coding the 
material properties and the interaction properties for each region. The first step is to find, for 
every nodal point, the neighboring nodal point. Those skilled in the art will appreciate that 
the operation of finding the neighboring nodal point corresponds to building the Voronoi 
diagram of the mesh. Its dual, the Delaunay triangulation, represents the best possible 
tetrahedral finite element for a given nodal configuration. The Voronoi diagram is 
constructed by a dilation approach. Under that approach, each nodal point in the discrete 
volume is dilated. Such dilation achieves two purposes. First, it is tested when one dilated 
point contacts another, so that neighboring points can be identified. Second, every voxel can 
be associated with a point of the mesh. 

Once every point x* has been associated with a neighboring point x,-, the two points are 
considered to be attached by a spring having spring constant kjj, where / and m identify the 

materials. The spring constant is defined by the material interaction properties of the 
connected points; those material interaction properties are predefined by the user in 
accordance with known properties of the materials. If the connected points belong to the 
same region, the spring constant reduces to fc/j and is derived from the elastic properties of 

the material in the region. If the connected points belong to different regions, the spring 
constant is derived from the average interaction force between the materials at the boundary. 

In theory, the interaction must be defined between any two adjacent regions. In 
practice, however, it is an acceptable approximation to define the interaction only between 
major anatomical components in the image and to leave the rest as arbitrary constants. In 
such an approximation, the error introduced is not significant compared with other errors 
introduced in the assumptions set forth above. 
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Spring constants can be assigned automatically, particularly if the region of interest 
includes tissues or structures whose approximate size and image intensity are known a priori, 
e.g., bone. Segmented image regions matching the a priori expectations are assigned to the 
relatively rigid elastic constants for bone. Soft tissues and growing or shrinking lesions are 
5 assigned relatively soft elastic constants. 

Once the mesh has been set up, the next image in the sequence is input at step 309, 
and the energy between the two successive images in the sequence is minimized at step 311. 
The problem of minimizing the energy U can be split into two separate problems: 
minimizing the energy associated with rigid motion and minimizing that associated with 
10 deformable motion. While both energies use the same energy function, they rely on different 
strategies. 

The rigid motion estimation relies on the fact that the contribution of rigid motion to 
the mesh deformation energy (AX r KAX)/2 is very close to zero. The segmentation and the a 
priori knowledge of the anatomy indicate which points belong to a rigid body. If such points 
15 are selected for every individual rigid region, the rigid motion energy minimization is 
accomplished by finding, for each rigid region Ri, the rigid motion rotation R/ and the 
translation T, that minimize that region's own energy: 
AAT r/gW =min Al t/^= ^(AX-mm^ [/.(A*,)) 

Vlerigid 

where AX/ = Rf X< + T/X/ and Ax, is the optimum displacement matrix for the points that 

20 . belong to the rigid region R(. That minimization problem has only six degrees of freedom for 
each rigid region: three in the rotation matrix and three in the translation matrix. Therefore, 
the twelve components (nine rotational and three translational) can be found via a six- 
dimensional steepest-descent technique if the difference between any two images in the 
sequence is small enough. 



15 



WO 03/025837 PCT/US02/29005 

Once the rigid motion parameters have been found, the deformational motion is 
estimated through minimization of the total system energy U. That minimization cannot be 
simplified as much as the minimization of the rigid energy, and without further 
considerations, the number of degrees of freedom in a 3D deformable object is three times the 
number of node points in the entire mesh. The nature of the problem allows the use of a 
simple gradient descent technique for each node in the mesh. From the potential and kinetic 
energies, the Lagrangian (or kinetic potential, defined in physics as the kinetic energy minus 
the potential energy) of the system can be used to derive the Euler-Lagrange equations for 
every node of the system where the driving local force is just the gradient of the energy field. 
For every node in the mesh, the local energy is given by 

(«(*.(*, + A*) -g n M)) 2 4- (/?(|Vg„(x, + Ax)| -|Vg B+1 (x f .)|)) 2 + 
r(V 2 g w + Ax) + V 2 g w+1 (x,)) 2 + ±* Yflfibi " * ~ 

where G m represents a neighborhood in the Voronoi diagram. 

Thus, for every node, there is a problem in three degrees of freedom whose 
minimization is performed using a simple gradient descent technique that iteratively reduces 
the local node energy. The local node gradient descent equation is 
x f (w + 1) = x.(«) - vAU Mn)tn) (Ax) 

where the gradient of the mesh energy is analytically computable, the gradient of the field 
energy is numerically estimated from the image at two different resolutions, x(w+l) is the 
next node position, and v is a weighting factor for the gradient contribution. 

At every step in the minimization, the process for each node takes into account the 
neighboring nodes' former displacement. The process is repeated until the total energy 
reaches a local minimum, which for small deformations is close to or equal to the global 
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minimum. The displacement vector thus found represents the estimated motion at the node 
points. 

Once the minimization process just described yields the sampled displacement field 
AX, that displacement field is used to estimate the dense motion field needed to track the 
5 segmentation from one image in the sequence to the next (step 313). The dense motion is 
estimated by weighting the contribution of every neighbor mode in the mesh. A constant 
velocity model is assumed, and the estimated velocity of a voxel x at a time / is v(x, 0 = 
Ax(t)/At. The dense motion field is estimated by 

v(*,o=~7r 2j i r 

A* VAx y eG ffl <*,) I*-*/ 1 

10 where 

-l-i 

c(x) = 



li ,m is the spring constant or stiffiiess between the materials / and m associated with the voxels 
x and x;, M is the time interval between successive images in the sequence, |x - Xy| is the 
simple Euclidean distance between the voxels, and the interpolation is performed using the 
1 5 neighbor nodes of the closest node to the voxel x. That interpolation weights the contribution 
of every neighbor node by its material property k$\ thus, the estimated voxel motion is 

similar for every homogeneous region, even at the boundary of that region. 

Then, at step 315, the next image in the sequence is filled with the segmentation data. 
That means that the regions determined in one image are carried over into the next image. To 
20 do so, the velocity is estimated for every voxel in that next image. That is accomplished by a 
reverse mapping of the estimated motion, which is given by 

1 v-i 

v(x,/ + A0 = — 2>(x ; ,0 

ft V[x y +v( X/ ,/)] 6 5(x) 
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where H is the number of points that fall into the same voxel space S(x) in the next image. 
That mapping does not fill ail the space at time f+Af, but a simple interpolation between 
mapped neighbor voxels can be used to fill out that space. Once the velocity is estimated for 
every voxel in the next image, the segmentation of that image is simply 

5 Z(x, t + At) = L(x - v(x, t + Af)At, t) 

where L(x,t) and L(x,t+At) are the segmentation labels at the voxel x for the times t and t+At. 

At step 317, the segmentation thus developed is adjusted through relaxation labeling, 
such as that done at steps 211 and 215, and fine adjustments are made to the mesh nodes in 
the image. Then, the next image is input at step 309, unless it is determined at step 319 that 

10 the last image in the sequence has been segmented, in which case the operation ends at step 
321. 

The operations described above can be implemented in a system such as that shown in 
the block diagram of Fig. 4. System 400 includes an input device 402 for input of the image 
data, the database of material properties, and the like. The information input through the 

15 input device 402 is received in the workstation 404, which has a storage device 406 such as a 
hard drive, a processing unit 408 for performing the processing disclosed above to provide 
the 4D data, and a graphics rendering engine 410 for preparing the 4D data for viewing, e.g., 
by surface rendering. An output device 412 can include a monitor for viewing the images 
rendered by the rendering engine 410, a further storage device such as a video recorder for 

20 recording the images, or both. Illustrative examples of the workstation 304 and the graphics 
rendering engine 410 are a Silicon Graphics Indigo workstation and an Mx Explorer 3D 
graphics engine. 

Shape and topology of the identified biomarkers can be quantified by any suitable 
techniques known in analytical geometry. The preferred method for quantifying shape and 
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topology is with the morphological and topological formulas as defined by the references 
cited above. 

The data are then analyzed over time as the individual is scanned at later intervals. 
There are two types of presentations of the time trends that are preferred. In one class, 
5 successive measurements are overlaid in rapid sequence so as to form a movie. In the 
complementary representation, a trend plot is drawn giving the higher order measures as a 
function of time. For example, the mean and standard deviation (or range) of a quantitative 
assessment can be plotted for a specific local area, as a function of time. 

The accuracy of those measurements and their sensitivity to subtle changes in small 

10 substructures are highly dependent on the resolution of the imaging system. Unfortunately, 
most CT, MRI, and ultrasound systems have poor resolution in the out-of-plane, or "z" axis. 
While the in-plane resolution of those systems can commonly resolve objects that are just 
under one millimeter in separation, the out-of-plane (slice thickness) is commonly set at 
1.5mm or even greater. For assessing subtle changes and small defects using higher order 

15 structural measurements, it is desirable to have better than one millimeter resolution in all 
three orthogonal axes. That can be accomplished by fusion of a high resolution scan in the 
orthogonal, or out-of-plane direction, to create a high resolution voxel data set (Pefla, J.-T., 
Totterman, S.M.S., Parker, KJ. "MRI Isotropic Resolution Reconstruction from Two 
Orthogonal Scans," SPIE Medical Imaging, 2001, hereby incorporated by reference in its 

20 entirety into the present disclosure). In addition to the assessment of subtle defects in 
structures, that high-resolution voxel data set enables more accurate measurement of 
structures that are thin, curved, or tortuous. 

In following the response of a person or animal to therapy, or to monitor the 
progression of disease, it is desirable to accurately and precisely monitor the trends in 

25 biomarkers over time. That is difficult to do in conventional practice since repeated scans 

19 



WO 03/025837 PCT/US02/29005 

must be reviewed independently and the biomarkers of interest must be traced or measured 
manually or semi-manually with each time interval representing a new and tedious process 
for repeating the measurements. It is highly advantageous to take a 4D approach, such as was 
defined in the above-cited patent to Parker et al, where a biomarker is identified with 

5 statistical reasoning, and the biomarker is tracked from scan to scan over time. That is, the 
initial segmentation of the biomarker of interest is passed on to the data sets from scans taken 
at later intervals. A search is done to track the biomarker boundaries from one scan to the 
next. The accuracy and precision and reproducibility of that approach is superior to that of 
performing manual or semi-manual measurements on images with no automatic tracking or 

10 passing of boundary information from one scan interval to subsequent scans. 

While a preferred embodiment of the invention has been set forth above, those skilled 
in the art who have reviewed the present disclosure will readily appreciate that other 
embodiments can be realized within the scope of the present invention. For example, any 
suitable imaging technology can be used. Therefore, the present invention should be 

1 5 construed as limited only by the appended claims. 



i 
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We claim: 

1. A method for assessing a cancerous tissue in a patient, the method comprising: 

(a) taking at least one three-dimensional image of a region of interest of the patient, 

the region of interest comprising the cancerous tissue; 
5 (b) identifying, in the at least one three-dimensional image, at least one biomarker of 

the cancerous tissue; 

(c) deriving at least one quantitative measurement of the at least one biomarker; and 

(d) storing an identification of the at least one biomarker and the at least one 
quantitative measurement in a storage medium. 

10 2. The method of claim 1, wherein step (d) comprises storing the at least one three- 

dimensional image in the storage medium. 

3. The method of claim 1, wherein step (b) comprises statistical segmentation of the at 
least one three-dimensional image to identify the at least one biomarker. 

4. The method of claim 1, wherein the at least one three-dimensional image comprises 
15 a plurality of three-dimensional images of the region of interest taken over time. 

5. The method of claim 4, wherein step (b) comprises statistical segmentation of a . 
three-dimensional image selected from the plurality of three-dimensional images to identify 
the at least one biomarker. 

6. The method of claim 5, wherein step (b) further comprises motion tracking and 
20 estimation to identify the at least one biomarker in the plurality of three-dimensional images 

in accordance with the at least one biomarker identified in the selected three-dimensional 
image. 

7. The method of claim 6, wherein the plurality of three-dimensional images and the 
at least one biomarker identified in the plurality of three-dimensional images are used to form 
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a model of the region of interest and the at least one biomarker in three dimensions of space 
and one dimension of time. 

8. The method of claim 7, wherein the biomarker is tracked over time in the model. 

9. The method of claim 1, wherein a resolution in all three dimensions of the at least 
one three-dimensional image is finer than 1 mm. 

10. The method of claim 1, wherein the at least one biomarker is selected from the 
group consisting of: 

• tumor surface area; 

• tumor compactness (surface-to-volume ratio); 

• tumor surface curvature; 

• tumor surface roughness; 

• necrotic core volume; 

• necrotic core compactness; 

• necrotic core shape; 

• viable periphery volume; 

• volume of tumor vasculature; 

• change in tumor vasculature over time; 

• tumor shape, as defined through spherical harmonic analysis; ' 

• morphological surface characteristics; 

• lesion characteristics; 

• tumor characteristics; 

• tumor peripheral characteristics; 

• tumor core characteristics; 

• bone metastases characteristics; 

• ascites characteristics; 
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• pleural fluid characteristics; 

• vessel structure characteristics; 

• jieovasculature characteristics; 

• polyp characteristics; 

• nodule characteristics; 

• angiogenisis characteristics; 

• tumor length; 

• tumor width; and 

• tumor 3d volume. 

11. The method of claim 1, wherein the quantitative measure is at least one of tumor 
shape, tumor surface morphology, tumor surface curvature and tumor surface roughness. 

12. The method of claim 1, wherein step (a) is performed through magnetic resonance 
imaging. 

13. A system for assessing a cancerous tissue in a patient, the system comprising: 

(a) an input device for receiving at least one three-dimensional image of a region of 
interest of the patient, the region of interest comprising the cancerous tissue; 

(b) a processor, in communication with the input device, for receiving the at least one 
three-dimensional image of the region of interest, identifying, in the at least one three- 
dimensional image, at least one biomarker of the cancerous tissue and deriving at least one 
quantitative measurement of the at least one biomarker; 

(c) storage, in communication with the processor, for storing an identification of the at 
least one biomarker and the at least one quantitative measurement; and 

(d) an output device for displaying the at least one three-dimensional image, the 
identification of the at least one biomarker and the at least one quantitative measurement. 
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14. The system of claim 13, wherein the storage also stores the at least one three- 
dimensional image. 

15. The system of claim 13, wherein the processor identifies the at least one 
biomarker through statistical segmentation of the at least one three-dimensional image. 

16. The system of claim 13, wherein the at least one three-dimensional image 
comprises a plurality of three-dimensional images of the region of interest taken over time. 

17. The system of claim 15, wherein the processor identifies the at least one 
biomarkers through statistical segmentation of a three-dimensional image selected from the 
plurality of three-dimensional images. 

18. The system of claim 17, wherein the processor uses motion tracking and 
estimation to identify the at least one biomarker in the plurality of three-dimensional images 
in accordance with the at least one biomarker identified in the selected three-dimensional 
image. 

19. The system of claim 18, wherein the plurality of three-dimensional images and the 
at least one biomarker identified in the plurality of three-dimensional images are used to form 
a model of the region of interest and the at least one biomarker in three dimensions of space 
and one dimension of time. 

20. The system of claim 13, wherein a resolution in all three dimensions of the at least 
one three-dimensional image is finer than 1 mm. 

21. The system of claim 13, wherein the at least one biomarker is selected from the 
group consisting of: 

• tumor surface area; 

• tumor compactness (surface-to-volume ratio); 

• tumor surface curvature; 

• tumor surface roughness; 
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• necrotic core volume; 

• necrotic core compactness; 

• necrotic core shape; 

• viable periphery volume; 

• volume of tumor vasculature; 

• change in tumor vasculature over time; 

• tumor shape, as defined through spherical harmonic analysis; 

• morphological surface characteristics; 

• lesion characteristics; 

• tumor characteristics; 

• tumor peripheral characteristics; 

• tumor core characteristics; 

• bone metastases characteristics; 

• ascites characteristics; 

• pleural fluid characteristics; 

• vessel structure characteristics; 

• neovasculature characteristics; 

• polyp characteristics; 

• nodule characteristics; 

• angiogenisis characteristics; 

• tumor length; 

• tumor width; and 

• tumor 3d volume. 

22. The system of claim 13, wherein the quantitative measure is at least one of tumor 

shape, tumor surface morphology, tumor surface curvature and tumor surface roughness. 

25 



WO 03/025837 



Fig, 1 



1/3 



PCT/US02/29005 



Take 3D image of region of interest 



Identify at least one biomarker in 
image and take quantitative measure 



Take multiple 3D images overtime 



Identify biomarker in the images taken 
over time and take quantitative measure 



I 



Develop 4D (space and time) model of 
region of interest with biomarker 



Track development of biomarker over 
time 



( o *f 



6 -£ 



"1 



fft?~ H INPUT { - 



==3h. 



) STORAGE | 



I 



PROCES SING ^t 08 
I 



GRAPHICS 
RENDERING 



I 



I 



OUTPUT H 4!2 



LWDRKmnpN^J 



WO 03/025837 



2/3 



PCT/US02/29005 



TAKE IMAGES "U ^) 

INPUT IMAGE | , 
(FIRST IN SEQUENCE) 'V'J' 



MEAN AND VARIANCE 
ESTIMATION VIA 
THE LINE SEGMENT 
ESTIMATOR 



CONNECTIVITY ESTIMATION 
VIA MEAN AND VARIANCE 
PAIR COMPARISON TEST 



I 



J 



REGION SPLITTING | 



FIG. 2 



GIBBS RANOOM FIELD BASED 
SEGMENTATION RELAXATION 



I 



STATISTICAL REGION 
MERGING VIA LOCAL 
BHATTACHARYYA METRIC 



I 



SEGMENTATION RELAXATION! 



WO 03/025837 



3/3 



PCT/US02/29005 



si. 



7 

LU 



to 

LU 



x ce 

LU LU 



LU LU 













as 








LU »-H 
























LU U_ 








yj 








LU J— 








h-« LL 




s: cd 




•— • 

















CO 




H- 
















5 


3 


LJ 


LU 


-< 




_J 




Q- 


□ 


tO 




h- * 


1— 


a 


o 




s: 


LU 




a 


LU 


LD 


CO 




-2L 




LU 

a 


s 






LU 






2: 




H-t 




h- 


LU 


1^3 


s 


LU 




HE 










LU 


LD 




2: 




>— 1 








3 



LU 

PC 



L3L 

ld 



rsiLO 

SB 

32 LU 



>- 
CD 



1 — l 

LU K— 
ID -aC 
<C 2£ 

1-4 CD 
LU 
CD Z 
LU 



lu a 
Ql- 

LUC3 
to 

LULU 

to 

LU CD 



3i£ 

CL. 

CD LD 

CL. >-* 
tO 




LD 



l-ltO LU 
LU LD 
LU D£ -<C 
LD at 3E 

2Z L_J 

LU LU 

E?E£ 

QC LU 

Har 

U.2LD 
LU 

Y- CXI 
ZD to 

a. lu 
to 



a 

CO I — 
LU -< 
HH CK 
1 — LU 

UJ tO 
CL. LU 
CD to I— r 
B2 LU >— 
CL. Q£ C£ 
LU 

UJ p- CL. 
Z3 t-J CD 

to zd a* 
to cv: a- 

h- tO 









NTS BY RANDDHL 
GRADIENT PDIN1 
INTS. RANODHU 
ERIOR POINTS'' 






»— 1 -CD H- 

a pr cl. z 

CL LD f-H 




1— « >- 
LU DC Q* UJ 

S LD CD CD 




DEFINE 
SELECT] 
ANO 601 
SELECT 



St 
32 

LU 

*g{ 

Lu S 
a a 

CD LU 
— I Li. 



1^1 

u. 



INTERNATIONAL SEARCH REPORT 



International application No. 
PCT/USQ2/290O5 



A. CLASSIFICATION OF SUBJECT MATTER 

IPC(7) : G06K9/00 

USCL : 382/131 
According to International Patent Classification (IPO or to both national classification and IPC 



B. FIELDS SEARCHED 



Minimum documentation searched (classification system followed by classification symbols) 
U.S. : 382/131 



Documentation searched other than minimum documentation to the extent that such documents are included in the fields searched 



Electronic data base consulted during the international search (name of data base and, where practicable, search terms used) 



DOCUMENTS CONSIDERED TO BE RELEVANT 



Category * 



Citation of document, with indication, where appropriate, of the relevant passages 



Relevant to claim No. 



X 
Y 
Y 
Y 
Y 
A 



US 6,1 12, 1 12 A (GILHUILS et al.) 29 August 2000 col. 3 lines 15-19; col. 4 lines 21-42 
and 33-35; col. 5 lines 1-10, col. 6 lines 60-63; col. 9 lines 38-53; and col. 10 lines 10- 
20. 

US 6,277,074 Bl (CHATURVEDI et al.) 21 August 2001 col. 4 lines 16-30 

US 5,785,654 A (IINUMA et al.) 28 July 1998 col. 8 lines 1-36 

US 6,368,331 Bl (FRONT et al.) 09 April 2002 col. 3 lines 44-55 

US 5,068,788 A (GOODBNOUGH et al.) 26 November 1991 col. 2 lines 62-65, col . 3 
lines 20-40, col. 6 lines 55-62 



1-5,10-12,13-17,21, 
and 22 



6-9, 18,19, and 20 
6-8,18, and 19 

6-8,18, and 19 

9 and 20 

1-22 



| | Further documents are listed in the continuation of Box C. See patent family annex. 



* Special categories of cited documents: 

*A" document defining the general itate of the art which Is not considered to be 
of particular relevance 

"E" earlier application or patent published on or after the international filing date 

«L» document which may throw doubts on priority ciaim(s) or which is cited to 
establish the publication date of another citation or other special reason (as 
specified) 

"O" document referring to an oral disclosure, use, exhibition or other means 

*P* document published prior to the international filing date but later than the 
priority date claimed 



"X" 



later document published after the International filing date or priority 
date and not in conflict with the application but cited to understand the 
principle or theory underlying the invention 

document of particular relevance; the claimed invention cannot be 
considered novel or cannot be considered to involve an Inventive step 
when the document Is taken alone 

document of particular relevance; the claimed invention cannot be 
considered to involve an inventive step when the document Is 
combined with one or more other such documents, such combination 
being obvious to a person skilled in the art 

document member of the same patent family 




Date of the actual completion of the international search 
15 November 2002 (15. 1 1 .2002) 



Name and mailing address of the ISA/US 

Commissioner of Patents and Trademarks 
Box PCT 

Washington, D.C. 20231 

Facsimile No. (703)305-3230 



Authorized offw 
Amelia Au 
Telephone No.l 103-'. 



Form PCT/ISA/210 (second sheet) (July 1998) 



